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Three  underground  cylindrical  structures  subjected  to 
an  impinging  time-dependent  plane  wave  have  been  analyzed  in 
this  study.  The  stress  wave  is  assumed  to  result  from  an 
underground  explosion.  The  structural  configurations  treated 
are  an  unlined  cavity,  a reinforced  concrete-steel  liner  com- 
posite cylinder  and  a cellular  concrete-steel  liner  composite 
cylinder.  The  medium  surrounding  the  structures  is  a soil-like 
material  (NTS  tuff)  and  the  impinging  wave  amplitude  (approxi- 
mately 0.5  kbar)  is  such  that  the  tuff  experiences  significant 
nonlinear  deformation.  Also,  both  composite  structures  respond 
in  the  nonlinear  regime. 


The  highly  nonlinear  nature  of  this  problem  requires  the 
use  of  numerical  simulation  techniques.  The  SWIS  finite  element 
code  (Frazier  and  Petersen,  1974;  Sweet,  ^ al . , 1976),  was 
utilized  because  of  its  ability  to  treat  both  the  structure 
and  the  surrounding  medium  in  the  nonlinear  response  regime. 

This  code  is  described  in  detail  in  Appendix  A.  Material  models 
utilized  for  this  study  include  the  Soil  Cap  Model  (DiMaggio 
and  Sandler,  1971)  for  the  NTS  tuff,  a variable  modulus  model 
for  the  reinforced  concrete  (Nelson,  e;t  al . , 1971)  , a von  Mises 
plasticity  model  for  the  steel  liner  and  a P-a  crush-up  model 
(Herrmann,  1969;  Carroll  and  Holt,  1972;  Good,  1972;  Riney, 
et  al . , 1972)  for  the  cellular  concrete. 


The  structural  configurations  considered  in  this  study 
are  meant  to  simulate  experiments  fielded  in  an  underground 
explosion  environment.  They  represent  three  different  design 
philosophies;  the  unlined  cavity  utilizes  the  inherent  strength 


of  the  medium,  the  reinforced  concrete  structure  strengthens 
the  cavity,  and  the  cellular  concrete  decouples  the  liner 
behavior  from  that  of  the  surrounding  medium.  The  following 
sections  discuss  the  results  of  the  SWIS  simulations. 
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CHAPTER  II 

PROBLEM  CONFIGURATION 

The  configurations  analyzed  in  this  study  are  the  response 
of  buried  cylindrical  structures  subjected  to  an  explosion- 
induced  stress  wave.  Effects  due  to  the  earth's  free  surface 
are  ignored.  The  structure  is  assumed  to  be  far  enough  removed 
from  the  explosion  so  that  a plane  wave  assumption  is  valid. 

Also,  end  effects  of  the  cylinders  are  not  treated  and  a plane 
strain  cross-sectional  configuration  is  utilized. 

Various  nomenclature,  including  the  configuration  of 
the  structure  and  surrounding  medium,  are  depicted  in  Figure  1. 

As  can  be  seen,  a plane  of  symmetry  allows  one-half  of  the 
problem  to  be  considered.  The  plane  wave  is  assumed  to  travel 
along  the  horizontal  direction  from  the  left  towards  the  right 
of  Figure  1. 

A finite  element  simulation  of  the  buried  structures 
requires  definitions  of  the  finite  element  mesh,  material  models, 
and  boundary  conditions.  The  following  siabsections  describe 
the  information  appropriate  for  the  simulations  of  this  study. 

2.1  FINITE  ELEMENT  REPRESENTATION 

The  finite  element  representation  of  the  dashed  region 
of  Figure  1 can  be  seen  in  Figure  2.  The  finite  element  mesh 
of  the  tuff  medium  as  well  as  that  of  the  two  composite  struc- 
tures are  depicted.  For  the  unlined  cavity,  all  elements 
surrounding  the  5-foot-diameter  cylindrical  opening  are  quad- 
rilateral. (This  unlined  cavity  mesh  is  similar  to  the  one 
appearing  in  Figure  B.l.) 

A study  similar  to  that  represented  by  Figure  1 has 
been  completed  for  a linear  structure-medium  configuration 
(Sweet,  1975B) . A finite  element  mesh  identical  to  that  of 
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(a)  Structure  configuration  and  finite  element  grid  location 


(b)  Angular  definitions. 


Figure  1.  Nomenclature  and  configuration  of  the  unlined 
cavity  and  composite  structures. 
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Figure  2 was  utilized  for  this  calculation.  The  excellent 
correlation  with  the  theoretical  solution  obtained  for  this 
study  serves  as  a validation  of  the  capability  of  SWIS  for 
treating  structure-medium  interaction  problems  of  this  type. 
This  linear  solution  is  also  summarized  in  Appendix  A. 

2.2  NTS  TUFF  MATERIAL  BEHAVIOR 

The  material  model  used  for  the  NTS  tuff  is  the  Soil 
Cap  Model  (DiMaggio  and  Sandler,  1971) . A detailed  description 
of  this  constitutive  model  can  be  found  in  Appendix  A.  The 
highly  nonlinear  nature  of  this  model  allows  the  observed 
compaction  behavior  of  this  soil-like  medium  to  be  simulated. 
The  material  parameters  used  in  this  model  are  given  below 
(Baron,  1975;  Sandler,  1975A;  Sandler,  1975B) : 

A = 0.5  kbar 

B = 0.52  kbar  ^ 

C = 0.44  kbar 

D = 1.8  kbar  ^ 

R = 3.0 

W = 0.015 

X = 0.003  kbar. 
o 


In  addition,  the  tuff  has  an  initial  density  of  2.0  gm/cm'^  and 
''' -Stic  bulk  and  shear  moduli  of  100  kbars  and  38  kbars,  respec- 
tively. The  tensile  strength  of  the  material  is  assumed  to  be 
zero.  The  procedure  used  in  SWIS  is  to  test  the  principal 
stresses  for  tensile  failure  before  the  cap  model  is  utilized. 
The  state  of  stress  is  first  assumed  to  be  elastic.  If  any 
principal  stress  is  tensile  it  is  set  to  zero.  The  resulting 
state  of  stress  is  then  modified  by  the  cap  model  if  necessary. 
This  method  is  superior  to  the  technique  of  using  (= 


) 

( 

il 
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as  the  criterion  for  tensile  failure  since  can  be  compressive 
and  yet  an  individual  principal  stress  can  be  tensile. 

The  uniaxial  strain  behavior  of  this  material  utilizing 
the  cibove  parametric  values  is  depicted  in  Figure  3.  Also 
appearing  in  this  figure  is  a comparison  with  the  uniaxial 
strain  behavior  for  the  material  properties  of  Baron,  1975. 

The  stress-strain  behaviors  for  hydrostatic  and  triaxial  states 
of  stress  show  similar  correlations  with  their  corresponding 
results  of  this  reference. 

In  the  actual  buried  structural  configuration,  of  course, 
the  medium  is  subjected  to  dji  situ  tectonic  and  gravitational 
stresses.  A dynamic  calculation  with  a static  prestressed 
gravitational  load  and  a stress-free  unlined  cavity  is  described 
in  Appendix  B.  The  inclusion  of  the  ^ situ  stresses  affected 
the  calculated  maximum  cavity  closure  by  approximately  15  per- 
cent for  this  configuration.  Since  this  difference  is  small 
compared  to  the  unknowns  involved  in  nonlinear  constitutive 
modeling  and  since  the  actual  ^ situ  stress  level  is  unknown, 
the  iji  situ  stresses  are  ignored  in  the  results  reported  in 
this  study. 


2.3 


CELLULAR  CONCRETE 


The  cellular  concrete  composite  structure  analyzed  in 
this  study  (designed  by  Lindberg,  1975)  requires  a cellular 
concrete  with  a 1500  psi  (0.1034  kbar)  crush  strength.  The 
actual  field  test  environment  is  the  dynamic  loading  of  a 
cellular  concrete  with  approximately  250-day  strength  properties. 
Existing  laboratory  material  properties  data  consist  of  static, 
28-day  values.  It  has  been  previously  noted  that  the  behavior 
of  cellular  materials  is  highly  strain-rate  dependent  (Butcher, 
1971;  Hoff,  1975).  It  is  also  true,  as  with  other  concrete 
materials,  that  the  behavior  of  cellular  concrete  depends  upon 
its  age.  Due  to  a lack  of  extensive  strain-rate  data,  a 
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STRESS  (kbar) 


Figure  3.  Uniaxial  strain  behavior  of  NTS  tuff  used 
in  calculations  (bold  line) . Behavior  of 
Baron  (1975)  (light  line)  is  also  given. 


rate-dependent  model  such  as  that  given  by  Butcher  (1971)  can- 
not be  formulated.  Rate-independent  properties  appropriate 
for  the  expected  strain  rate  will  be  utilized  instead.  Using 
Equations  (8)  and  (9)  as  given  by  Hoff  (1975) , yields  a static, 
28-day  strength  of  900  psi  (0.062  kbar)  if  a dynamic,  250-day 
strength  of  1500  psi  is  desired.  The  P-a  crush-up  model  described 
in  Appendix  A will  be  used  to  represent  the  irreversible  crush- 
up  behavior. 

It  will  be  assumed  that  the  "strength"  refers  to  the 
value  obtained  in  a uniaxial  strain  environment.  Using  avail- 
able experimental  data  (Hoff,  1971;  Hoff,  1972;  Hoff,  1975), 
the  stress-volumetric  strain  response  up  to  a maximum  strain 
of  0.47  is  given  in  Figure  4.  These  data  are  the  result  of 
punch  tests  and  thus  are  not  exactly  uniaxial  strain;  however, 
the  differences  will  be  ignored.  The  P-a  model  requires  P(v) 
data  where  P is  the  pressure  and  v is  the  specific  volume. 

Uniaxial  strain  data  may  be  used  to  derive  P(v)  if  the  shear 
behavior  of  the  porous  material  is  known.  It  will  be  assumed 
that  the  deviatoric  stress  behavior  of  the  cellular  concrete 
is  described  by  a simple  von  Mises  yield  stress,  Y^,  and  the 

associated  flow  rule.  The  shear  modulus,  y , is  assumed  to 

o 

be  constant.  Thus,  for  a uniaxial  strain  environment  (e^^  0) 

with  principal  stresses  (positive  in  comparison)  > O2  = cr^, 
the  linear  and  plastic  deviatoric  stress  behavior  is  given  by 


(A) 

linear:  *^i  ” *^2  ~ 

(la) 

(b) 

plastic:  ~ ^2 

= Y 

0 

(lb) 

either 

case  the  pressure  is 

defined  as 

P = 

°1  ^^2 

(2) 

3 
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Assuming  that  and  are  known,  Equation  (A-10)  of  Appendix 
A and  Equations  (1)  and  (2)  may  be  used  with  the  experimental 
data  of  Figure  4 to  derive  a P(v)  relationship.  Using  avail- 
able data,  it  is  assumed  that  the  porous  bulk  modulus  and  the 
porous  density  are  given  by: 

K = 9 kbar 
o 

(3a) 

= 0.854  gm/cm 

and  from  Lee  (1973),  the  solid  density  is  given  by: 

p =2.25  gm/cm^  . (3b) 

“o 

For  concrete  without  aggregate  it  may  be  assumed  that  the 
solid  bulk  modulus  becomes  (Read,  1971) ; 

K = 200  kbar  . (3c) 

s 

Furthermore,  it  is  assumed  that  the  Poisson's  ratio  of  the 
cellular  material  is  0.38  yielding  a shear  modulus  (using 
Equation  (3a) 

= 2.4  kbar  . (3d) 

Finally,  the  yield  strength  is  assumed  to  be 

Y = 0.024  kbar  . (3e) 

o 

The  data  of  Equation  (3)  along  with  Figure  4 are  sufficient  to 
prescribe  the  28-day  behavior  of  the  cellular  concrete  up  to 
a maximum  stress,  of  4500  psi  (0.31  kbar).  For  consistency 

a crush  pressure  of  0.88  kbar  is  assumed.  Using  the  P-a 
expressions,  data  of  Figure  4,  data  of  Equation  (3),  and 
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assuming  = 0.88  kbar,  the  resulting  P(v)  behavior  appears  \ 

in  Figure  5.  Since  it  is  well-known  that  strength  and  modulus  j 

of  concrete  increase  at  approximately  the  same  ratio,  the  250-  j 

day,  dynamic  values  are  determined  by  scaling  the  28-day  prop-  j 

erties  by  the  stress  ratio  1500/900  e 5/3.  This  P(v)  curve  j 

also  is  given  in  Figure  5.  The  resulting  a(P)  relationship  is  | 

given  in  Figure  6.  ] 

i 

The  effect  of  water  saturation  on  the  behavior  of  porous  I 

materials  is  extremely  important  (Riney,  et  al . , 1972) . It  j 

can  be  expected  that  the  material  behavior  of  Figures  5 and  6 ' 

would  be  significantly  altered  if  the  cellular  concrete  were  I 

saturated.  j 

2.4  REINFORCED  CONCRETE  MATERIAL  BEHAVIOR 

The  reinforced  concrete  material  behavior  is  modeled 
using  a variable  modulus  model,  (Nelson,  et  al . , 1971) . Its 
formulation  is  described  in  Appendix  A.  Parametric  values 
for  the  material  constants  are  given  by: 

hydrostatic  material 
properties 

shear  material 

f 

parameters  [ 

i 
5 

3 ' 

In  addition,  its  initial  density  is  2.415  gm/cm  . The  above  j 

properties  result  in  an  8000  psi  compressive  strength.  This 

strength  has  been  increased  from  the  28-day  static  value  of 
5500  psi  to  account  for  dynamic  and  aging  effects. 
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= -100 
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Figure  5.  P(v)  relationship  for  28-day  (static)  and 
250-day  (dynamic)  cellular  concrete. 
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Figure  6.  a(P)  relationship  for  28-day  (static)  and 
250-day  (dynamic)  cellular  concrete. 
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STEEL  MATERIAL  PROPERTIES 


The  material  model  used  for  the  inner  steel  liner  of 
the  composite  structures  is  a von  Mises  yield  model  utilizing 
the  associated  flow  rule.  This  plasticity  model  is  described 
in  detail  in  Appendix  A.  The  following  material  properties 
for  the  steel  are  used: 

Density  = 7.84  gm/cm 
Bulk  Modulus  = 1,590  kbar 
Shear  Modulus  = 775  kbar 

2.6  STRESS  BOUNDARY  CONDITION 

The  plane  wave  stress  boundary  condition  of  Figure  1 is 
applied  at  the  left  boundary  of  the  mesh  of  Figure  2.  The 
resulting  stress  and  velocity  time  histories  across  from  the 
structure  at  the  bottom  of  the  mesh  as  predicted  by  SWIS  are 
depicted  in  Figure  7.  Also  appearing  in  this  figure  are  the 
behaviors  suggested  by  Baron,  1975.  The  steepness  of  the 
SWIS  results  are  due  to  the  fact  that  a more  refined  mesh  was 
used  by  SWIS.  This  same  stress  boundary  condition  was  used 
to  perform  the  iji  situ  stress  study  of  Appendix  B. 


Yield  Stress  = 2.48  kbar 
(36  ksi) 

Young's  Modulus  = 2,000  kbar 
Poisson's  Ratio  = 0.29 
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Figure  7.  Free-field  stress  and  velocity  transients  at  structure  horizontal 


CHAPTER  III 
CALCULATIONAL  RESULTS 


The  dynamic  behaviors  of  the  ualined  cavity  and  the 
two  composite  structures  are  discussed  in  this  section. 

Data  typically  measured  in  a field  experiment  are  presented. 

The  finite  element  meshes  depicted  in  Figure  2 were  utilized 
and  the  impinging  plane  wave  description  is  depicted  in  Figure  7. 
Velocity  histories  at  several  points  in  the  tuff  and  the  struc- 
ture are  presented.  Also,  the  hoop  strain  behavior  in  the 
composite  structures  are  also  included.  An  important  factor 
that  can  be  used  to  assess  the  behavior  of  these  structures  is 
the  cavity  closure,  which  is  inherently  related  to  the  stability 
of  these  underground  openings.  Cavity  closure  time  histories 
are  presented  for  the  three  structural  configurations. 


3.1  UNLINED  CAVITY 

Three  velocity  time  histories  on  the  horizontal  axis 
of  the  unlined  cavity  calculation  can  be  seen  in  Figure  8. 

The  nomenclature  for  the  angular  orientation  is  depicted  in 
Figure  1.  The  horizontal  velocity  on  the  vertical  axis  appears 
in  Figure  9.  The  obvious  reflections  from  the  cavity  can  be 
seen  in  both  these  figures.  Dynamic  cavity  closure  is  depicted 
in  Figure  10.  As  can  be  seen,  closure  occurs  on  both  the 
horizontal  and  vertical  axes. 

3.2  REINFORCED  CONCRETE-STEEL  LINER  COMPOSITE  STRUCTURE 

Velocity  time  histories  in  the  tuff  surrounding  the 
reinforced  concrete  structure  appear  in  Figures  11  and  12. 

The  reflections  from  this  structure  are  seen  to  be  less  pro- 
nounced compared  to  those  of  the  unlined  cavity.  Horizontal 
velocity  behavior  at  three  locations  on  the  inside  of  the 
steel  liner  appears  in  Figure  13.  Increased  high  frequency 
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Figure  9.  Horizontal  velocity  waveform  in  the  tuff  on  the 
vertical  axis,  approximately  12  feet  from  the 
unlined  cavity  (i.e.,  90°  or  270°,  12'). 
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Figure  10.  Diametrical  behavior  of  the  unlined  cavity. 
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Figure  12.  Horizontal  velocity  waveform  in  the  tuff  on  the 
vertical  axis,  approximately  12  feet  from  the 
reinforced  concrete  structure  (i.e.,  90®  or 
270®,  12’)  . 
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content  is  noticeable  for  the  steel  behavior  compared  to  the 
free-field  input. 


The  0.5  kbar  stress  wave  results  in  nonlinear  behavior 
in  both  the  reinforced  concrete  and  steel  liner,  as  well  as 
the  surrounding  tuff.  Hoop  strain  in  the  structural  components 
resulting  from  the  impinging  plane  is  depicted  in  Figures  14  and 
15.  As  can  be  seen,  the  focusing  of  energy  around  the  struc- 
ture is  particularly  apparent  in  the  steel  liner.  The  cavity 
closure  appearing  in  Figure  16  exhibits  closure  on  the  hori- 
zontal axis  and  opening  on  the  vertical  axis. 

3.3  CELLULAR  CONCRETE-STEEL  LINER  COMPOSITE  STRUCTURE 

Because  of  the  low  compression  strength  of  the  cellu- 
lar concrete,  the  velocity  behavior  in  the  tuff  surrounding 
this  structure,  which  appears  in  Figures  17  and  18,  is  quite 
similar  to  the  unlined  cavity  results  of  Figures  8 and  9. 

Also,  the  reduced  strength  of  the  concrete  has  modified  the 
steel  liner  velocity  histories,  presented  in  Figure  19,  com- 
pared to  those  of  the  reinforced  structure.  Figure  13. 

Strains  in  the  cellular  concrete  and  steel  liner  are 
given  in  Figures  20  and  21.  Comparing  Figures  15  and  21, 
quite  dissimilar  behavior  in  the  steel  liner  is  noted  for 
the  two  composite  structures.  The  cavity  closure  depicted 
in  Figure  22  is  similar  to  the  reinforced  structure  in  that 
closure  on  the  horizontal  axis  and  opening  on  the  vertical 
axis  occurs. 
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Figure  14. 
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Hoop  strain  waveforms  in  the  center  of  the  reinforced 
concrete  at  the  45°  (or  315°)  and  135°  (or  225°) 
locations. 
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Figure  15.  Hoop  strain  waveforms  at  the  steel  inner  surface 
of  the  reinforced  concrete  structure  at  the  45° 
(or  315°)  and  135°  (or  225°)  locations. 
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Figure  16.  Diametrical  behavior  of  the  reinforced  concrete 
structure.  Arrow  in  the  inset  figure  is  the 
direction  of  the  plane  wave  and  lies  along  the 
horizontal  axis. 
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Figure  18.  Horizontal  velocity  waveform  in  the  tuff  on  the 
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Figure  20.  Hoop  strain  waveforms  in  the  center  of  the 
cellular  concrete  at  the  45°  (or  315°)  and 
135°  (or  225°)  locations. 
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Figure  21.  Hoop  strain  waveforms  at  the  steel  inner  surface 
of  the  cellular  concrete  structure  at  the  45° 

(or  315°)  and  135°  (or  225°)  locations. 
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Figure  22.  Diametrical  behavior  of  the  cellular  concrete 

structure.  Arrow  in  inset  figure  is  the  direction 
of  the  plane  waves  and  lies  along  the  horizontal 
axis . 
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CHAPTER  IV 

SUMMARY  AND  DISCUSSION 

Results  presented  in  the  previous  section  serve  to 
illustrate  the  difference  in  the  three  designs  analyzed  in 
this  study.  The  unlined  cavity  does  not  attempt  to  increase 
the  strength  inherent  in  the  surrounding  tuff;  while  the 
reinforced  concrete  approach  relies  on  the  strength  of  the 
composite  structure  for  integrity.  However,  as  illustrated 
in  Figure  15,  the  steel  lining  experiences  a significant 
strain  environment  in  the  reinforced  design  concept.  The 
use  of  cellular  concrete  reduces  the  maximum  strain  in  the 
steel  liner  by  approximately  a factor  of  three. 

Maximum  cavity  closure  also  varies  significantly 
among  the  three  cavity  configurations.  These  values  are 
approximately  7.5  cm,  3.8  cm  and  2 cm  for  the  unlined  cavity, 
reinforced  composite  and  cellular  composite,  respectively. 

The  cellular  concrete  approach  is  obviously  more  effective 
since  it  decouples  the  behavior  of  the  interior  of  the  struc- 
ture from  the  behavior  of  the  surrounding  tuff  medium.  This 
important  advantage  of  the  cellular  concrete  design  can  become 
even  more  significant  for  media  with  compressive  strengths 
greater  than  tuff.  For  this  type  of  medium,  the  reinforced 
design  must  also  be  strengthened  to  be  effective.  The  dis- 
advantage of  this  concept  is  that  it  may  become  prohibitively 
impractical  to  strengthen  concrete  much  beyond  8000  psi.  It 
should  be  cautioned  that  the  effectiveness  of  the  cellular 
concrete  can  be  reduced  if  it  becomes  saturated  with  water. 
This  is  an  important  fact  considering  the  wet  environment  of 
underground  structures. 
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SWIS  FINITE  ELEMENT  COMPUTER  CODE  - NONLINEAR  VERSION 


The  SWIS  Finite  element  computer  code  (Frazier  and 
Petersen,  1974  and  Sweet,  e^  al . , 1976)  considers  either 
static  or  dynamic  problems  in  one-,  two-,  and  three-dimensional 
geometries  with  equal  efficiency.  In  addition,  both  nonlinear 
and  linear  behavior  may  be  analyzed.  The  dynamic  calculations 
may  be  performed  using  either  the  explicit  or  implicit  time 
integration  schemes.  Also,  for  explicit  dynamic  calculations, 
SWIS  has  the  capability  of  varying  the  time  step  within  the 
computational  mesh.  This  capability  has  been  quite  useful  for 
structure-medium  calculations  (Sweet,  1976A  and  Sweet,  1976B) . 
The  present  discussion  primarily  concerns  the  nonlinear  version 
of  the  SWIS  code.  However,  comparisons  with  known  linear 
theoretical  solutions  are  also  included  for  information  pur- 
poses . 


A.l  INTRODUCTION 

SWIS  has  been  formulated  in  a manner  that  is  similar 
to  Lagrangian  finite  difference  techniques  in  several  respects. 
The  spatial  discretization  is  performed  at  each  computational 
step  to  preclude  the  formation  of  a stiffness  matrix.  Also, 
stress  is  computed  from  strain  in  an  incremental  fashion 
explicitly  in  each  element  at  each  computational  step.  The 
resulting  computational  procedure,  therefore,  has  the  advantage 
of  treating  irregular  geometries  as  well  as  the  advantages  of 
minimal  storage  requirements  and  the  ease  of  incorporating 
nonlinear  constitutive  relationships. 

The  equilibrium  equations  are  derived  using  conventional 
finite  element  techniques.  Interpolation  within  skewed  elements 
is  accomplished  in  an  isoparametric  manner.  This  interpolation 
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scheme  is  then  combined  with  the  virtual  work  expression  to 
obtain  the  discrete  governing  equations.  The  mass  matrix  is 
diagonalized  for  computational  efficiency.  Element  integrals 
are  performed  using  a single-point  quadrature.  Bending  within 
an  element  is  modeled  by  utilizing  the  exact  solution  to  simple 
bending  (both  linear  and  nonlinear  bending  solutions  are  applied 
when  appropriate) . This  modified  single-point  quadrature  tech- 
nique has  been  demonstrated  to  be  valid  for  several  problem 
configurations  (see  Section  A. 5).  a summary  of  the  features 
of  SWIS  are  as  follows: 

• Finite  element  formulation 

• One-,  two-,  or  three-dimensional  spatial  con- 
figurations 

• Isoparametric  elements  with  superior  bending 
characceri sties 

• Lagrangian  large  displacement  formulation 

• Choice  of  several  nonlinear  material  models 

• Dynamic  calculations  using  explicit  central 
difference  time  integration 

• Equilibrium  calculations  using  the  conjugate 
gradient  method  (static  or  dyneimic  implicit) 

• The  ability  in  a dyneunic  calculation  to  inte- 
grate a portion  of  the  finite  element  mesh 
(usually  the  structure)  using  a smaller  time 
step  than  in  the  surrounding  elements 

• Both  regular  and  arbitrary  element-node  num- 
bering 

• Restart  and  graphics  capabilities. 

A. 2 EXPLICIT  DYNAMIC  CALCULATIONS 

SWIS  employs  a time-centered  method  for  integrating  the 
equilibrium  equations  for  explicit  calculations.  Given  the 
state  of  stress  at  time,  t,  the  calculation  procedure  proceeds 
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(2)  The  velocity,  displacement  and  nodal  position 
vectors  are  updated  as  follows: 


u(t  + = u(t  - + dt  u(t) 

u(t  + dt)  = u(t)  + dt  u (t  + 

X (t  + dt)  = x(t)  + dt  u(t  + ^) 


(3)  Strain  rate  at  time  t + ^ is  determined  from 
the  velocity  vector. 

(4)  Stresses  at  time  t + dt  are  determined  from 
stresses  at  time,  t,  and  the  strain  rate  at 

time  t ^ using  nonlinear  constitutive  rela- 
tions. The  stress  tensor  is  also  modified  to 
account  for  material  rotations. 

(5)  Time  is  incremented  and  steps  (1)  - (4)  are 
repeated. 

A. 3 STATIC  AND  IMPLICIT  DYNAMIC  CALCULATIONS 

The  incremental  discrete  equilibrium  equations  for  a 
dyneimic  configuration  may  be  written  in  the  following  form 

mAvi”  = + Ar"  (A-1) 

n n“l 

where  Ag  = g - g , m is  the  mass  matrix,  vi  is  the 

acceleration  vector , F is  the  applied  force  vector,  R is 
the  restoring  force  vector,  and  the  superscript  n refers  to 
time  t”.  Time  integration  is  accomplished  by  using  either  the 
midpoint  constant  acceleration  algorithm  (Powell,  1972)  or  the 
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Newmark  6-method  (Nickell,  1973)  . In  either  case,  the  accel- 


eration increment  at  t*'  can  be  expressed  as  a function  of  the 


displacement  increment  and  the  motion  at  time  t”~^,  i.e.. 


A A^n  . n 1 n— 1 

Au  = a Au  + b 


(A-2) 


where  a is  a known  constant  and  b is  a function  of  the  motion 


ri"!  n 

at  time  t . A first  order  approximation  for  Ar"  is  given  by 


AR 


,n  _ /iR\ 

\9u/ 


n-1 


Au 


n 


(A-3) 


Combining  Equations  (A-1)  - (A-3)  yields 


n 


A Au  = B 


(A-4) 


where  A is  the  effective  stiffness  matrix  and  B is  the 
effective  restoring  force.  These  terms  are  defined  by 

n-1 


A = m a 


-(If) 


B = Af”  - m b"“^ 


Equation  (A-4)  is  solved  in  SWIS  using  the  conjugate  gradient 
method  (Marcus,  1960). 


n-1 


Since  the  stiffness  term  A is  evaluated  at  time  t* 
the  above  approach  yields  an  unbalanced  nodal  force  at  the  end 
of  each  time  step.  The  procedure  used  in  SWIS  to  obtain  the 
correct  solution  is  to  use  the  solution  of  Equation  (A-4)  as 
the  first  approximation,  ^^^Au”.  The  unbalanced  forces  are 
then  used  to  obtain  the  next  iteration  and  this  procedure  is 


repeated  until  convergence  is  noted  (Bathe,  et  al.  1974), 
This  iterative  procedure  for  the  k^^  iteration  is  summarized 


by  the  following  equation: 
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ii 


J . 

■ 

i 


1 


A(<Wiu" 


= Af"  + r" 


- R + m(u 


(k-1)  ..n. 
u ) 


where  equals  r”  In  order  to  make  this  iteration 

procedure  more  efficient,  the  secant  modulus  may  be  used  to 
define  the  matrix,  A.  Nonlinear  static  calculations  are  per- 
formed by  eliminating  the  inertial  terms  and  interpreting  t” 
as  load  step  n. 


A. 4.  MATERIAL  MODELS 


A. 4.1  Von  Mises  Plasticity  Model 

The  von  Mises  plasticity  model  is  a constitutive  rela- 
tionship that  models  a material  as  an  elastic-perfectly  plastic 
continuum.  The  deviatoric  stress  tensor,  ^ , is  related  to 
the  total  stress  tensor  by 


'kk 


A yield  function  is  defined  as 


Y 

/3 


where  Y is  the  one-dimensional  yield  stress.  The  deviatoric 
behavior  is  assumed  to  be  elastic  when 


F < 0 


and  perfectly  plastic  when 
F = 0. 
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According  to  the  associated  flow  rule,  the  strain  tensor. 


is  assumed  to  be  composed  of  elastic  and  plastic  com- 


ponents , 


e . . = ef . + e? . 
1]  3-1  3-1 


and  the  plastic  strain  rate  tensor  is  related  to  the  yield 
function  by 


9F 


e? . = X » 
ij  3 a.. 


3-1 


where  X is  a scalar  multiplier.  The  elastic  strain  rate 
tensor  is  assumed  to  be  related  to  the  state  of  stress  by 
Hooke's  law.  This  relationship  becomes 


¥) 


where  y is  the  shear  modulus  (assumed  to  be  constant) . As- 
siaming  a plastic  state  of  stress  (F  = 0)  yields  the  following 
expressions  for  X and  S, 


3-1 


X = 


/3  e.  . S.  . 
iJ il 


S.  . = 2y 
3-1 


kk 


i: 

3 


/T  X 


2Y 


’ij) 


For  completeness,  the  hydrostatic  state  of  stress  is  defined 
as  a function  of  the  state  of  strain.  In  nonlinear  analyses. 


the  pressure,  p (=  - » is  usually  represented  as  a poly- 


nomial function  of  (P/P^  “ 1)»  i.e.. 


" ‘ Si  (t  ■ 
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where  p is  the  density  and  the  density  ratio  is  related  to 
the  strain  by 

^ = exp  (-Ekk’  • 

The  propagation  of  a plane  wave  in  an  elastic-perf ectly  plastic 
material  has  been  simulated  in  SWIS.  The  result  of  this  calcu- 
lation is  depicted  in  Figure  A.l  where  the  correlation  is  seen 
to  be  excellent. 

A. 4. 2 Soil  Cap  Constitutive  Model 

The  soil  cap  model  (DiMaggio  and  Sandler,  1971;  Read  and 
Frazier,  1975;  and  Sandler  and  Rubin,  1976)  is  a plasticity  model 
whose  plastic  strain  rate  vector  is  normal  to  the  yield  sur- 
face in  stress  spaoe.  The  yield  surface  (depicted  in  Figure  A. 2) 
is  composed  of  a fixed  failure  envelope  that  encompasses  a 

movable  hardening  cap. 

Defining  and  J'  as; 

"^1  " ‘^kk 


the  yield  function  is  given  by 

= /J^  - [A  - C exp  (B  J^)]  , if  L < , 

= - L)^  + R^  J2  - (X  - L)^  , if  L > Ji 

where  X and  L are  defined  as 

[(1  - ]/D  * X„ 

Z for  £ ^ 0 
0 for  £ < 0 


X = £r 
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Figure  A.l.  Comparison  of  SWIS  and  theoretical  solutions  for  an  elastic-plastic 
uniaxial  strain  calculation. 


and  a satisfies 

I 

Jl  + R [A  - C exp  (-BJI)]  = X 

The  terms  A,  B,  C,  D,  R,  W and  Xq  are  material  constants. 

The  n\americal  algorithms  which  represent  the  above 
expressions  are  available  in  SWIS.  The  models  of  both  Read 
and  Frazier,  1975  and  Sandler  and  Rubin,  1976  are  available. 
Typical  stress-strain  behavior  of  this  model  can  be  seen  in 
Figures  A. 3 and  A. 4. 

An  example  SWIS  calculation  using  the  cap  model  is 
summarized  in  Figures  A. 5 and  A. 6.  Here,  the  propagation  of 
a plane  wave  through  a two  layered  soil  has  been  simulated 
(Read  and  Frazier,  1975) . As  can  be  seen,  the  correlation 
with  theory  is  excellent. 


i 

I 

f 

r, 

f 


A. 4. 3 Variable  Moduli  Model  , 

This  incremental  stress-strain  constitutive  model  assumes  ' 
that  the  material  behavior  can  be  separated  into  volumetric  i 
and  deviatoric  parts  with  no  explicit  yield  behavior  (Nelson,  i 
et  al . , 1971)  . This  separation  automatically  precludes 

dilatancy  in  the  material.  Both  the  bulk  (volumetric)  modulus  s 
and  shear  (deviatoric)  modulus  depend  upon  the  stress  and/or  | 
strain  invariants  and  the  previous  loading  history.  ] 

The  bulk  modulus  is  assumed  to  obey  the  following  rela-  ^ 
tionship;  j 
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Axial 


Figure  A. 3.  Typical  uniaxial  strain  behavior  of 
the  soil  Cap  model. 


Stress  Difference  (kbar) 


Figure  A. 5.  The  uniaxial  strain  response  of  the  clay  and  shale  used  in  two-layered 

example  calculation.  The  materials  are  simulated  using  the  soil  cap  model 


Stress  (ksi)  Particle  Velocity  (in/msec) 


0 5 10  15  20  25  30  35  40  45 


Distance  (feet) 

(a)  Particle  velocity  versus  distance 


Distance  (feet) 

(b)  Stress  versus  distance 


Figure  A. 6.  Comparison  of  theoretical  and  computational  results 
for  one-dimensional  two-layered  calculations. 
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K = 


2 , 

K-  + K,e  + K-e  when  p = p and  p > 0 
0 12 

Ko  P ' Pmax  P ' Pmax  i < <> 


where  Pj^^x  maximum  previous  pressure,  e is  the  mean 

compressive  strain  and  are  material  parameters. 

The  shear  modulus  behavior  is  given  by: 


G = 


^0  ^l'^  ^1  when  > 0 

Gq  + when  £ 0 


where  is  the  second  deviatoric  stress  invariant,  p is  the 
pressure,  and  (Gq,  Y-j^)  are  material  parameters. 

The  above  material  model  yields  a behavior  that  is  charac- 
terized by  hardening  in  shear  with  increasing  pressure  and 
softening  with  increasing  shear  stress  > 0 and  y^  ^ 0)* 

A . 4 . 4 P-g  Constitutive  Model 

The  hydrostatic  behavior  of  crushable  materials  may  be 
represented  using  the  "P-g"  model  (Herrmann,  1969;  Riney,  et  al . , 
1972;  Good,  1972  and  Carroll  and  Holt,  1972) . This  approach  for 
representing  the  crushup  of  porous  materials  assumes  that  the 
equation-of-state  of  the  solid  (i.e.,  the  non-porous  material) 
is  known.  If  the  distension  ratio,  a,  is  defined  as  the  ratio 
of  the  solid  and  porous  densities,  i.e.. 


and  the  equation-of-state  of  the  solid  is  given  by 


then  the  hydrostatic  behavior  of  the  porous  material  is  as- 
sumed to  be  given  by 

P = i f (ap)  . (A-5) 

The  P-a  formulation  is  completed  by  specifying  a functional 
relationship  between  a and  P;  a relationship  that  is  usually 
experimentally  determined.  Writing  this  expression  as 

a = a(P)  , {A-6) 


the  composite  pressure,  P,  can  be  determined  as  a function  of 
the  composite  density,  p,  by  eliminating  a from  Equations 
(A-5)  and  (A-6) . The  term  a may  be  related  to  the  porosity, 
by 


A typical  stress-strain  curve  for  crushable  materials 
is  depicted  in  Figure  A. 7.  As  can  be  seen,  it  is  characterized 
by  an  initial  nearly-linear  behavior  followed  by  a highly  non- 
linear character  as  the  material  crushes.  These  data  can  be 
utilized  in  order  to  formulate  a P-a  model  if  the  behavior  of 
the  solid  material  is  known.  Assuming  that  the  low  kilobar 
stress  regime  is  of  interest,  the  solid  equation-of-state  becomes 

^ = ■'s  -)  • 


Thus,  using  Equations  (A-5)  and  (A-7) , the  composite  pressure  is 
given  by 


or,  in  terms  of  the  specific  volume. 


and  the  initial 
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Stress 


[ 


Figure  A. 7.  Typical  stress-strain  curves  for  crushable 
materials. 


solid  density  p ; 

®0 


Solving  the  above  expression  for  a yields 


(A-8) 


(A- 9) 


p V K 


Therefore,  if  P (v)  is  known,  the  a(P)  relationship  follows 
from  Equation  (A-9) . For  example,  the  volumetric  strain  of 
Figure  A. 7 defined  as 


^ = 1 _ V , ^ ^ ^ ^ ^ ^ 

^0  ^0  ° 

Therefore,  v is  related  to  ^V/V.  by: 


(A-10) 


Equations  (A-9)  and  (A-10)  along  with  experimental  data  yield 
a tabular  expression  for  a(P).  Two  additional  parameters  that 
describe  the  crushup  behavior  are  the  elastic  limit,  P , and 
the  fully  crushed  pressure,  P . These  terms  along  with  the 
unloading  behavior  are  depicted  in  Figure  A. 8.  Upon  unloading, 
complete  recovery  of  pore  space  is  assumed  when  the  pressure 
is  below  the  elastic  limit.  In  this  regime,  rather  than  use 
the  experimental  a(P)  data,  it  is  assumed  that  a(P)  is  of 
the  following  form: 


a(P)  = aP  + bP  + a, 


(A-11) 


1 


1 


Figure  A. 8. 


Schematic  of  P-a  model  for  the  distension  ratio, 
a(P).  Note  that  above  the  elastic  limit,  release 
from  a crush  state  occurs  along  a prescribed 
(reversible)  release  path  which  could  result  in 
some  void  recovery  (02  > 0.2)  . Release  in  the 
elastic  regime  is  confined  to  the  loading  curve. 
After  exceeding  the  crush  pressure,  Pq,  the 
material  loads/unloads  according  to  the  bulk 
modulus  of  the  solid. 


a 


61 


is  defined  as 


where 


The  parameters  a and  b are  determined  from  the  initial 
composite  bulk  modulus/  and  the  values  of  P and  v at 

the  elastic  limit,  and  v^.  The  bulk  modulus  K is  defined 

as 


K = p 


dP 
0 dp 


Therefore,  using  the  relation 


P = 0 


and  Equation  (A-8)  yields 


Noting  that 


results  in  the  following  expression  for  b: 


The  parameter  a is  determined  by  evaluating  Equation  (A-11) 

at  P'  and  becomes 
e 
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a 


(A-14) 


!• 


I 


a 


- b P - Otn 
e 0 


where 


follows  from  Equation  (A-9) 


(A-15) 


When  the  pressure  is  between  the  elastic  limit  and  the 
crush  pressure,  the  tabular  expression  for  a(P)  is  used  with 
Equation  (A-8)  to  yield  P(v)  or  a(v).  Partial  pore  recovery 
occurs  when  the  pressure  unloads  from  this  regime.  If  a*  is 
the  minimum  value  of  a ever  obtained  and  a is  the  value  of 
a using  the  tabulated  a(v)  relation,  then  the  value  of  a 
during  unloading  is  given  by  (Good,  1972): 

a*-  a 

a = ^ (a*  - a)  + a . (A-16) 

“e 

Thus,  during  unloading,  the  amount  of  distension  that  is  re- 
covered is  assumed  to  be  linearly  proportional  to  the  excursion 
from  the  elastic  limit. 

Negative  pressures  are  not  allowed  by  the  following  tech- 
nique. If  a negative  pressure  is  calculated,  a new  value  of  a 
is  determined  by  setting  P of  Equation  (A-8)  to  zero.  This 
yields 


0 


or 


a = 


(A-17) 


I 


'i 


kL 
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In  sununary,  the  material  properties  required  for  the 
implementation  of  this  constitutive  model  are: 


• tabular  relation 

• initial  density  and  bulk  modulus;  and 


solid  properties,  p and  K 

Sq  s 


& 


where  the  tabular  Pl^l  yields  the  values  of  P^,  P^  and  v^. 


A. 4. 5 Tensile  Failure  Model 

There  exists  a large  class  of  problems  where  the  geometry 
and/or  loading  configuration  results  in  a significant  region  of 
interest  where  the  tensile  behavior  is  important.  One  modeling 
procedure  that  has  been  utilized  to  simulate  tensile  failure  is 
to  simply  zero  all  calculated  stresses  which  exceed  the  tensile 
strength  (Sandler,  et  al. , 1972) . This  technique  is  similar  to 
approaches  which  alter  the  material  stiffness  properties  in  that 
it  does  not  offer  any  insight  concerning  the  inelastic  strains 
(e.g. , crack  porosity)  which  result  from  tensile  failure.  The 
tensile  failure  model  described  here  is  an  extension  to  the 
model  presented  by  Maenchen  and  Sack,  1964  and  Cherry,  Sweet 
and  Halda,  1973.  Briefly,  inelastic  strains  are  introduced  in 
order  to  zero  those  principal  stresses  which  exceed  the  material 
tensile  strength.  The  accximulation  of  these  strains  is  the 
tensile  failure-induced  porosity.  In  the  present  analysis, 
careful  attention  is  paid  to  the  effect  of  this  porosity  on 
the  determination  of  the  mean  stress  from  the  equation  of  state. 

The  initial  step  in  the  computational  procedure  is  to 
calculate  the  state  of  stress  utilizing  the  equation  of  state 
for  the  mean  stress  and  assuming  elastic  behavior  to  derive 
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the  deviatoric  stresses.  For  a two-dimensional  (x,y)  config- 
uration, these  stresses  are  given  by  the  deviatoric  components 

PT  Using  standard  techniques. 


S'  ) and  the  mean  stress 
xy 

the  principal  stresses  are  given  by 


(S  ; S', 
X y 


T'  = - (P'  + S'  + S') 
11  X y 


"^22'  "^33  ■ ^ 2 


S'  + S'+  V(S'-S')^  + (2S' 

X y-  1 X y xy 


Each  of  these  stresses  is  compared  to  the  tensile  strength  of 
the  mateiral  if  no  previous  tensile  failure  has  occurred  in  the 
computational  element.  A zero  tensile  strength  is  appropriate 
if  previous  failure  has  taken  place.  When  one  of  the  principal 
stresses  exceeds  the  tensile  strength,  an  increment  of  strain 
in  this  direction  is  introduced  to  alter  the  primed  state  of 
stress.  This  strain  increment  is  accumulated  for  each  cycle 
and  the  sum  for  all  three  principal  directions  is  interpreted 
as  the  tensile  failure-induced  porosity.  The  strain  increment 
results  in  the  principal  stress  increments  (AT^^,  , AT^.,,  AT,,) 
and  yields  the  following  final  state  of  stress; 


T = T ' + AT 

^11  ^11  ^11 


"^22  "^22  ^*^22 


(A-18) 


T = T ' + AT 

33  33  33 


Upon  fracturing,  the  stress  increments  must  be  chosen  such  that 
the  principal  stress  that  exceeded  the  tensile  limit  is  "relaxed" 
to  zero  and  also  the  mean  stress  resulting  from  the  above  expres- 
sions is  consistent  with  the  equation  of  state.  This  stress 
adjustment  is  assumed  to  occur  with  the  deviatoric  stress  com- 
ponent determined  using  linear  elasticity.  Thus: 
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where  AP  is  the  increment  in  mean  stress,  y is  the  shear 

modulus,  and  As..  is  the  inelastic  strain  increment  in  the 
th 

i direction.  The  sign  of  As..  has  been  chosen  so  that  a 


positive  increment  results  in  a positive  change  in  the  porosity 
content.  The  double  subscripts  do  not  denote  summation.  Using 
an  approach  similar  to  that  developed  for  multi-phase  material 
modeling  (Riney,  ^ al . , 1972  and  Sweet,  1972),  the  mean  stress 
is  given  by 


V 

P = ^ f (A-20) 

V m 

where  V is  the  specific  volume  of  the  computational  cell, 

is  that  portion  of  V occupied  by  the  matrix  material 

(volume  of  voids  equals  , and  f is  the  equation  of  state 

of  the  poreless  matrix  material.  In  general,  the  equation  of 

state  is  assumed  to  be  a function  of  V and  the  internal 

m 

energy  per  unit  mass,  E. 

As  previously  mentioned,  the  fracturing  process  intro- 
duces an  increment  of  porosity  in  the  calculational  cell.  Since 
the  volume  of  the  cell  is  not  altered,  this  results  in  a change 
in  the  matrix  content.  Thus,  the  AP  term  required  in  Equation 
(A-19)  and  calculated  from  Equation  (A-20)  is  given  by 

= If-  * H “ • 

m 

The  use  of  Equation  (A-21)  ensures  that,  at  least  in  an  incre- 
mental sense,  the  adjusted  mean  stress  satisfies  the  equation 
of  state.  The  matrix  volume  change  is  related  to  the  inelastic 
strain  increments  by 
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V 


(A-22) 


AV  = 
m 


Ae , 


11 


i=l 


and  the  energy  increment/  AE,  can  be  related  to  these  strain 
increments  and  the  state  of  stress  using  the  energy  equilibrium 
equation.  This  equilibrium  relationship  may  be  written 


E = - P V + V S. .e. . 

1]  ID 


where  e. . is  the  total  strain  rate  tensor  and  S. . are  the 
ID  ID 

deviatoric  stresses.  Remembering  that  the  fracturing  process 

does  not  change  V,  the  above  expression  yields  the  following 

relationship : 


i=l  i=l 


Tii  + P 


Ae.  . 
' 11 


(A-23) 


When  tensile  failure  occurs,  any  number  of  the  primed 
principal  stresses  may  exceed  the  tensile  limit.  The  case  when 
one  of  these  stresses  exceeds  this  limit  will  be  investigated 
first.  If  this  stress  is  the  strain  increment  Aej^j^  must 

be  found.  The  stress  increment  equals  and  thus 


Equations 

(A-18) 

and 

(A-19)  become 

■^11 

= T ^ 
^11 

- AP 

4 

3 

uACj^j^  E 0 

’*22 

= t' 
22 

- AP 

UAEjj 

(A-24) 

"^33 

= t' 
^33 

- AP 

UAEi,  . 

Equations  (A-21)  - (A-23)  are  used  to  derive  the  expression  for 
AP.  This  relationship  becomes 


AP 


V 


3P  V 3P 

3V  V-  3E 

m 0 


(T£i  + 


p') 


Ae 


11 


(A-25) 
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where  the  partial  derivatives  and  pressure  are  evaluated  at  the 
primed  stress  state.  Equation  (A-25)  may  be  written 


“ ' ■'eft  ''"ll 


(A-26) 

where  the  effective  bulk  modulus  and  is  defined  by 


K 


ef  f 


V - ^ (T"  + P") 

V.  3E  '^11  ^ ^ ’ 
m 0 


Using  Equation  (A-26) , Equation  (A-24)  becomes 

) ^^11 
I ^^11 
I ^^11 


*^11  = 

^il  - (^eff 

'’^22  ' 

*^22  ■ (^eff 

2 

3 

■^33  = 

*^33  " (^eff 

2 

3 

using 

Equation  (A- 

■24)  ; 

40^  = 

T 

^11 
^ 4 , 

^eff  3 ^ 

(A-27) 


(A-28) 


When  large  tensile  stresses  are  allowed  before  fracturing  occurs, 
the  direct  application  of  Equation  (A-27)  and  (A-28)  can  result 
in  numerical  instabilities  if  too  large  a stress  drop  is  allowed 
on  one  cycle.  This  situation  can  be  prevented  by  reducing  the 
increment  allowed  by  Equation  (A-28)  for  the  first  few  cycles 
after  fracturing  is  initiated. 


When  two  principal  stresses  exceed  the  tensile  limit. 


and  the  strain  increments  Ae 


say  ciiiu  ■*•22' 

introduced  to  "relax"  both 


and 


11 

to  zero. 


and  Ac 


are 


22 

The  analysis 


11  "22 

is  similar  to  that  previously  presented  for  one  strain  increment. 
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For  this  case,  AP  is  given  by 


= -V  ^ / 
\ 


^^22 


* (’'22  ^ ''')''"22] 


(A-29) 


As  can  be  seen,  its  form  does  not  lend  itself  to  the  definition 
of  an  effective  bulk  modulus.  The  stress  state  is  now  defined 
by 


Til  = f£i  - AP 


T = T " - AP 

22  22 


T33  = T33  - AP  + 


- |w  ( 

- |w  ( 


2iCii  - AE32 


2AC22  - ^^11 


(A-30) 


A^ll  ^^22 


Using  Equation  (A-29)  and  the  first  two  equations  of  (A-30) 
yields  expressions  for  Ae,,  and  Ae__.  They  are: 


where  a and  b follow  from  rewriting  Equation  (A-29)  as 


AP  = aAe^j^  + hLz22 


and  are  given  by 
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When  all  of  the  principal  stresses  exceed  the  tensile 
limit,  the  material  has  become  rubble  and  all  of  the  stresses 
must  be  zeroed.  For  this  case,  the  strain  increments  are  all 
assumed  to  equal 
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(A-32) 


For  a plane  stress  configuration,  the  above  procedure  is  modi- 
fied by  introducing  a crack  porosity  normal  to  the  plane  whose 
value  is  determined  by  constraining  the  normal  stress  to  be  zero. 

The  analysis  just  presented  is  also  valid  for  crack  clo- 
sure. Once  failure  has  occurred  and  nonzero  crack  porosity 
exists.  Equation  (A-28) , (A-31)  or  (A-32)  is  used  regardless 
of  the  sign  of  the  principal  stresses.  The  existence  of  com- 
pressive stresses  automatically  decreases  the  accvunulated  crack 
strain.  This  process  is  continued  until  the  strain  is  detected 
to  be  negative.  When  this  occurs,  the  strain  component  is  equat- 
ed to  zero  and  the  material  is  assumed  to  be  "healed,"  with  a 
future  zero  tensile  strength. 


A. 5 EXAMPLE  CALCULATIONS  PERFORMED  USING  THE  SWIS  CODE 

The  following  figures  summarize  several  calculations 
that  have  been  performed  using  the  SWIS  code.  As  can  be  seen, 
comparison  with  theory  is  excellent  for  all  excimples.  Figures 
A. 9 through  A. 11  demonstrate  that  the  bending  characteristics 
of  the  SWIS  elements  are  valid  for  both  linear  and  nonlinear 
behavior.  Figure  A. 12  illustrates  that  the  SWIS  formulation 
is  also  valid  for  continuum  wave  propagation. 
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Percent  Error  in  End  Deflection 


CANTILEVER  BEAM 

X,  in  equilibrium  under  applie 

end  load 

•x 

(Poisson's  ratio  = 1/4) 

V 

conventional 
finite  element 
(bilinear  displacements) 


,SWIS 


Figure  A. 9, 


Number  of  Elements 


A comparison  of  the  bending  scheme  used  in  the 
SWIS  code  with  that  inherent  in  conventional 
bilinear  finite  elements  for  computing  the  end 
deformation  of  a cantilever  beam. 


suppcjrted 


24"  X, 
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Iq  cos  6 
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(a)  Geometric  problem  configuration. 
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(C)  Radial  displacement  at  X3  > 0 and  8 
(b)  Longitudinal  displacement  at  X3  ■ l*/2  and  as  a function  of  time. 

0 ■ 0 as  a function  of  time. 


Figure  A. 10.  Comparison  between  theoretical  and  SWIS  results 
for  3-D  shell  problem. 
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(a)  Simply  supported  beam  and 
applied  load. 


(b)  Finite  element  idealization 
(twenty-four  elements)  for 
coarse  mesh.  Fine  mesh  has 
36  elements  along  length 
and  3 elements  through 
thickness. 
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— ~ — Linear  Solution,  SHIS  and  Theoretical 
• • • SUIS  Solution,  fine  nesh 

O O O 5WIS  Solution,  Coarse  Kesh 


Kaximun  Deflection 
Baron,  etal. , 1961 


O o 


0 2 


A > Static  Clastic  Deflection  of 
Beam  Subjected  to  p 

0 


"a  e 10  12  14  16 


Time  Ratio  (At/0.3  x 10'^) 


(c)  Comparison  of  SWIS  solutions  with  theoretical  solutions. 
Fine  Mesh. 


Figure  A. 11.  SWIS  simulation  of  the  dynamic  response  of  a 

simply  supported  beam  obeying  a von  Mises  yield 
criterion. 
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The  subcycling  capability  is  particularly  important 
I for  dynamic  explicit  structure-media  calculations.  As  a 

y 

i test  of  this  procedure,  a calculation  was  performed  for  the 

[ dynamic  response  of  a buried  cylinder  subjected  to  a 

plane  wave  (Sweet,  1976B) . The  finite  element  grid  for  this 
problem  can  be  seen  in  Figure  A. 13.  The  characteristic  time 
step  in  the  cylinder  is  approximately  8 times  smaller  than  that 
in  the  surrounding  medium.  Two  calculations  were  performed: 
one  without  subcycling  utilizing  a uniform  time  step  in  the 
grid  and  another  calculation  with  a time  step  in  the  sub- 
cycled region  being  one-eighth  that  in  the  surrounding  medium. 

The  excellent  comparison  for  these  two  calculations  can  be  seen 
in  Figure  A. 14. 

The  example  summarized  in  Figure  A. 15  considers  the 

dynamic  response  of  a structure-medium  interaction  configuration. 

The  variables  y , V and  Oq  not  defined  in  this  figure  are 

c pm  0 

the  shear  modulus  in  the  cylinder,  the  p-wave  velocity  in  the 
medium  and  the  hoop  stress  in  the  cylinder,  respectively. 
Parametric  values  of  variables  defined  in  Figure  A. 15  are  as 
follows: 

A = 5,  = 0.2;  pulse  width  and  cylinder  geometry 


n»  _ 

E 

0.4; 

Young ' s modulus  ratio 

c 

V = 

m 

0.25 

> Poisson  ratios 

V = 

c 

0.2 

) 

Pm  “ 

2.68 

gm/cm^  \ 

[ densities 

Pc  = 

2.32 

gm/cm^  ) 

Figure  A. 14.  Comparison  of  subcycled  (dashed  line)  and  uniform  time  step  (solid  line) 
results . 
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Figure  A. 15.  SWIS  simulation  of  the  dynamic  response  of  a buried 

thick-walled  cylinder  subjected  to  an  impinging  plane 


r 


Figure  A. 16  summarizes  the  three-dimensional  simulation 
of  a buried  fault.  The  configuration  is  a (prescribed)  propa- 
gating 100  cm  dislocation  and  thus  does  not  consider  the 
rupture  process. 
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APPENDIX  B 

THE  EFFECT  OF  ^ SITU  GRAVITATIONAL  STRESSES 
ON  THE  BEHAVIOR  OF  AN  UNLINED  CAVITY 

The  analyses  presented  in  this  report  assume  that  the 
explosion-induced  stress  wave  acts  in  a stress-free  medium. 

In  reality,  of  course,  ^ situ  overburden  as  well  as  tectonic 
stresses  exist  in  the  neighborhood  of  the  cavity.  A precise 
simulation  of  the  behavior  of  the  cavity  would  first  consider 
the  actual  fabrication  of  the  cavity  and  the  resulting  relief 
of  stresses  near  the  cavity  wall.  The  stress  wave  would  then 
be  treated  as  a stress  superimposed  on  this  existing  state 
of  stress.  Because  the  tuff  behaves  in  a nonlinear  fashion, 
the  response  of  the  cavity  can  be  expected  to  differ  somewhat 
depending  upon  whether  or  not  the  iji  situ  stresses  are  con- 
sidered . 

The  effect  of  ^ situ  gravitational  stresses  on  the 
response  of  an  unlined  cavity  is  discussed  in  this  appendix. 
This  simulation  is  accomplished  by  prestressing  the  external 
boundaries  of  a two-dimensional  finite  element  mesh  containing 
a cylindrical  cavity.  The  external  stress  consisted  of  an 
88  bar  vertical  stress  and  an  18  bar  horizontal  stress.  A 
stress-free  boundary  condition  was  maintained  for  the  cavity 
wall.  After  this  nonlinear  static  state  of  stress  was  imposed, 
a dynamic  stress  wave  boundary  condition  was  applied  at  the 
left  hand  mesh  boundary  (see  Figure  B.l)  to  simulate  an 
impinging  plane  wave.  It  should  be  noted  that  both  the  static 
and  explicit  dynamic  capabilities  of  the  nonlinear  version 
of  the  SWIS  code  were  utilized.  The  tuff  material  model  and 
the  dynamic  stress  wave  description  are  identical  to  those 
discussed  in  Section  2. 


Results  for  this  calculation  and  for  a calculation 
without  gravitational  effects  can  be  seen  in  Figures  B.2 


i 
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Plane  of  Symmetry 


Figure  B.l.  Finite  element  mesh  and  boundary  conditions  utilized  to  simulate  the  effect  of 
in  situ  overburden  stresses. 


through  B.5.  Figure  B.2  contains  a comparison  of  the  free 
field  (i.e.,  at  right  boundary  across  from  cavity)  stress  and 
velocity  for  these  two  calculations.  As  can  be  seen,  the 
peak  stresses  are  quite  similar  and  the  peak  values  of  velo- 
city differ  by  approximately  12  percent.  The  comparisons 
of  the  horizontal  components  of  velocity  at  the  inside  surface 
of  the  cavity  appearing  in  Figures  B.3  and  B.4  also  exhibit 
approximately  12  percent  differences.  Figure  B.5  contains  the 
comparison  of  the  horizontal  cavity  closure  for  the  calcula- 
tions with  and  without  gravitational  effects.  A fifteen  per- 
cent reduction  in  the  cavity  closure  when  gravity  is  included 
can  be  noted. 
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Velocity  (cm/ms) 


Time  (ms) 


Figure  B.3.  Horizontal  velocity  comparison  at  the  inside  of  the 
unlined  cavity  at  the  0°  (closest  to  left  boundary) 
location. 
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Velocity  (cm/ms) 


Time  (ms) 


Figure  B.4.  Horizontal  velocity  comparison  at  the  inside  surface 
of  the  unlined  cavity,  90“  location. 
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Diametrical  closure  (cm) 


Figure  B.5.  Comparison  of  the  cavity  closure  in  the  horizontal 
direction  for  the  calculations  with  and  without 
gravitational  effects. 
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